function [Tmat, wholemat, renewmat, fossilmat] = mainfun(x)


thetaj = x(1);
thetaB = x(2);
Tbar = x(3);
k_Tbar = x(4);
S_T0 = x(5);

global parms i kr Hr jr cr Sf kf Nf cf gf sigmaf T0
y0 = 1.0;     % output of goods y0 - choose units so this is 1.0
%
Sbar = 2126.0527;    % Total feasible resources
Res_0 = 1./0.065100642;       % Initial estimate of proved reserves
rho_0 = Sbar - Res_0;    % Initial ratio of Alpha2 to Alpha3
g_0 = 0.11068716;      % Initial value of per unit mining cost
n0 = 0.00832913;      % Initial level of mining investment
%
% % gN_0 = -0.081639946;    % Partial derivative of g with respect to N at t=0
% % Alpha3 = -gS_0.*Rho_0./gN_0;

alpha3 = 15; %-gS_0.*Rho_0./gN_0;
alpha2 = rho_0.*alpha3;
%
delta = 0.04;
%
c0 = 0.6619974; % Initial value of consumption for calculate lambda_0
k0 = 3.6071282734; % Initial value of K for the differential equation
%
A =  y0./k0;
%
% Marginal cost of backstop energy p = (Gamma1+H)^(-alpha)
% Alpha is the slope of the learning curve
%
alpha = 0.25;
%
% Given alpha, Gamma1 determines the initial cost of renewable energy. Here we set it
% to 4 times the initial cost of fossil fuel.
%
Gamma1 = (4*g_0).^(-1./alpha);
%
% % After some time t, marginal cost will decline to Gamma2 and remain there.
% % We assume this ultimate minimum marginal cost of renewables
% % is 20% of the initial cost p when H = 0:
%
Gamma2 =0.8*g_0;
%
Abar = A*(1-Gamma2)+ (1-delta);
% Note: We should have Abar > 0. It is different from Abar in continuous
% time model

% To make discretized model comparable to continuous one, beta is redefined
% so that the long term growth rate of the two models are the same.
% (beta*Abar)^(1/gamma)-1 = 4.07%, from which beta = 0.9673
% beta_new = 1/(1+beta), where beta = 0.05.
beta  =  0.9524;

% psi is the effect of learning by doing on H
% psi has to be between 0 and 1
% A smaller value of psi will allow a larger role for explicit investment
% in renewable technology as opposed to learning by doing

psi = 0.33;
psr = 1/psi;
psc = 1-psi;
%
% Q = population growth that is used to "scale" resource extraction
% (variables other than fossil fuel exploitation are in per capita terms)
% popgr is the exogenous population growth rate

Q0 = 1;
popgr = 0.01;
%
% % gamma is the coefficient of relative rsik aversion. If we wish
% % to calibrate to a particular initial consumption level, we can
% % allow gamma to vary.
gamma = 4.0;
gamr = 1/gamma;

gS_0 = 0.00015;  % Partial derivative of g with respect to S at t=0
alpha1 = gS_0.*Res_0.^2;
alpha0 = g_0 - alpha1./Res_0;

parms = [delta A Sbar alpha0 alpha1 alpha2 alpha3 Gamma1 alpha psi beta ...
    gamma Q0 popgr Abar];
% % Note that parms is a 1x14 vector



%%%%%%%%%%%%%%%  Regime 1: B>0, j>=0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ======================================
% At Tbar, the marginal cost of backstop technology becomes Gamma2. Hence
% we have (Gamma1+H)^(-alpha)=Gamma2
H_Tbar = Gamma2^(-1/alpha)-Gamma1;

% Calculate k_t, H_t and V_t for all t in the renewable regime.

% Initialize the solution matrix
tt = 500; % tt should be larger than # of years in renewable regime
ir = zeros(tt,1);
jr = zeros(tt,1);
cr = zeros(tt,1);
kr = zeros(tt,1);
Hr = zeros(tt,1);
etar = zeros(tt,1);
lambdar = zeros(tt,1);
PRr = zeros(tt,1); %% p_t+1/p_t




% Value of variables at Tbar
kr(1) = k_Tbar;
Hr(1) = H_Tbar;
ir(1) = ((beta*Abar)^gamr-1+delta)*k_Tbar;
cr(1) = (1-Gamma2)*A*k_Tbar-ir(1);
PRr(1) = 1/Abar;
etar(1) = 0;
lambdar(1) = cr(1)^(-gamma);
i = 2;
while true
    
    PRr(i) = 1/(A-(Gamma1+Hr(i-1))^(-alpha)*A+1-delta+(psi-thetaB)*jr(i-1)/((psc-thetaj)*kr(i-1)));
    cr(i) = (beta*cr(i-1)^(-gamma)/PRr(i))^(-gamr);
    lambdar(i) = cr(i)^(-gamma);
    Dbar = (psc-thetaj)*alpha*A*kr(i-1)*(Gamma1+Hr(i-1))^(-alpha-1)+A^(-psi)*kr(i-1)^(-psi)*jr(i-1)^psi;
    Cbar = Dbar*PRr(i);
    
    j = @(k) Cbar.^psr*A*k;
    inv = @(k) kr(i-1)-(1-delta)*k;
    H = @(k) Hr(i-1)-Cbar^(psr*psc)*A*k;
    
    if H(kr(i-1))>0 && inv(kr(i-1))>0
        kr(i) = fzero(@(k) (1-(Gamma1+H(k)).^(-alpha))*A.*k-inv(k)-j(k)-cr(i),kr(i-1));
        Hr(i) = H(kr(i));
        
        
        jr(i) = j(kr(i));
        ir(i) = inv(kr(i));
        etar(i) = beta*(etar(i-1) + lambdar(i-1)*alpha*(Gamma1+Hr(i-1))^(-alpha-1)*A*kr(i-1));
        i = i+1;
    else
        
        Hr(i) = 0;
        kr(i) = fzero(@(k) (1-Gamma1.^(-alpha))*A.*k-inv(k)-j(k)-cr(i),kr(i-1));
        jr(i) = j(kr(i));
        ir(i) = inv(kr(i));
        etar(i) = beta*(etar(i-1) + lambdar(i-1)*alpha*(Gamma1+Hr(i-1))^(-alpha-1)*A*kr(i-1));
        % Cut out the zero lines of solution matrix
        yearsr = i;
        Hr = Hr(1:yearsr);
        kr = kr(1:yearsr);
        jr = jr(1:yearsr);
        ir = ir(1:yearsr);
        cr = cr(1:yearsr);
        PRr = PRr(1:yearsr);
        etar = etar(1:yearsr);
        lambdar = lambdar(1:yearsr);
        break
    end
    
end




%%   %%%%%%%%%%%%%%%  Regime 2: R>0, n>=0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Key state variables: k and N
% Control variables: i and n
% Other state variables: Q and S
% Transition point T0




g = @(S,N) alpha0 + alpha1./(Sbar-alpha2./(alpha3+N)-S);
gdS = @(S,N) alpha1*(alpha3+N).^2./((Sbar-S).*(alpha3+N)-alpha2).^2;
gdN = @(S,N) -alpha1.*alpha2./((Sbar-S).*(alpha3+N)-alpha2).^2;


T0 = Tbar - yearsr;

% Set the initial results
invf = zeros(T0+1,1);
nf = zeros(T0+1,1);
kf = zeros(T0+1,1);
Nf = zeros(T0+1,1);
Sf = zeros(T0+1,1);
sigmaf = zeros(T0+1,1);
cf = zeros(T0+1,1);
Q = zeros(T0+1,1);
gf = zeros(T0+1,1);
PRf = zeros(T0+1,1);
nuf = zeros(T0+1,1);
lambdaf = zeros(T0+1,1);

kf(1) = kr(end);
cf(1) = cr(end);
lambdaf(1) = lambdar(end);
invf(1) = ir(end);
PRf(1) = PRr(end);

Q(1) = (1+popgr)^(T0);
gf(1) = Gamma1^(-alpha)-psi*etar(end)*(A*kf(1))^(psi-1)*jr(end)^(1-psi)/(cf(1)^(-gamma));
Sf(1) = S_T0;
Nf(1) = alpha2/(Sbar-Sf(1)-alpha1/(gf(1)-alpha0))-alpha3;


z_T0 = [kf(1),Sf(1),Nf(1),lambdaf(1),sigmaf(1),nuf(1)];
tspan0 = [T0, T0-1];
% %
% options = odeset('RelTol',1e-6,'AbsTol',1e-6,'events',@events0);
% [tf,zf,TE0,ZE0,IE0] = ode45(@foss,tspan0,z_T0,options);
%
options = odeset('RelTol',5e-14,'AbsTol',5e-14);
[tf,zf] = ode45(@foss,tspan0,z_T0,options);

% Extract components of z for graphing

kf(2)=zf(end,1);
Sf(2)=zf(end,2);
Nf(2)=zf(end,3);
lambdaf(2)=zf(end,4);
sigmaf(2)=zf(end,5);

Q(2) = (1+popgr)^(T0-1);
cf(2) = lambdaf(2)^(-gamr);
gf(2) = g(Sf(2),Nf(2));
invf(2) = inv(kf(2));
PRf(2) = (beta*lambdaf(1))/lambdaf(2);
nf(2) = Nf(1)-Nf(2);

for i = 3:T0+1
    
    Q(i) = (1+popgr)^(T0-i+1);
    cf(i) = (beta*((1-gf(i-1)+sigmaf(i-1)*Q(i-1)/(cf(i-1)^(-gamma)))*A+1-delta))^(-gamr)*cf(i-1);
    lambdaf(i) = cf(i)^(-gamma);
    
    sigmaf(i) = beta*(sigmaf(i-1)-gdS(Sf(i-1),Nf(i-1))*A*kf(i-1)*cf(i-1)^(-gamma));
        
    x0 = [kf(i-1),Nf(i-1)] ;
    options = optimset('Display','off','TolFun',1e-8,'TolX',1e-8);
    xopt = fsolve(@FOC_fossil,x0,options);
    kf(i) = xopt(1);
    Nf(i) = xopt(2);
    
    nf(i) = Nf(i-1)-Nf(i);
    invf(i) = kf(i-1)-(1-delta)*kf(i);
    Sf(i) = Sf(i-1)-Q(i)*A*kf(i);
    gf(i) = g(Sf(i),Nf(i));
    % PRf(i) =1/(1-A*kf(i-1)*gdN(Sf(i-1),Nf(i-1)));
    PRf(i) = (beta*lambdaf(i-1))/lambdaf(i);
    
end

P_Er = (Gamma1+Hr).^(-alpha)-psi.*etar.*(A.*kr).^(psi-1).*jr.^(1-psi)./lambdar;
P_Ef = gf - sigmaf.*Q./lambdaf;

disp(['Tbar = ',num2str(Tbar),';']);
disp(['k_Tbar = ', num2str(k_Tbar,18), ';']);
disp(['S_T0 = ',num2str(S_T0,12),';']);
disp(['T0 = ', num2str(T0),', k0 = ', num2str(kf(end),7), ',N0 = ',num2str(Nf(end)),', S0 = ',num2str(Sf(end))]);

Tmat = [Tbar, T0];
wholemat = [[P_Er;P_Ef],[PRr;PRf],[lambdar;lambdaf],[ir;invf],[cr;cf],[kr;kf]];
wholemat = reshape(wholemat(end:-1:1),Tbar+1,6);

renewmat = [etar jr Hr P_Er PRr lambdar ir cr kr];
renewmat = reshape(renewmat(end:-1:1),yearsr,9);
fossilmat = [sigmaf gf nf Nf Sf P_Ef PRf lambdaf invf cf kf];
fossilmat = reshape(fossilmat(end:-1:1),T0+1,11);
